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Abstract 

The ability of one- and two-equation turbulence 
models to predict unsteady separated flows over air- 
foils is evaluated. An implicit, factorized, upwind- 
biased numerical scheme is used for the integration of 
the compressible, Reynolds averaged Navier-Stokes 
equations. The turbulent eddy viscosity is obtained 
from the computed mean flowfield by integration 
of the turbulent field equations. The two-equation 
turbulence models are discretized in space with an 
upwind-biased, second order accurate total variation 
diminishing scheme. One and two-equation turbu- 
lence models are first tested for a separated airfoil 
flow at fixed angle of incidence. The same models 
are then applied to compute the unsteady flowfields 
about airfoils undergoing oscillatory motion at low 
subsonic Mach numbers. Experimental cases where 
the flow has been tripped at the leading edge and 
where natural transition was allowed to occur nat- 
urally are considered. The more recently developed 
field-equation turbulence models capture the physics 
of unsteady separated flow significantly better than 
the standard k — e and k — to models. However,- cer- 
tain differences in the hysteresis effects are obtained. 
For an untripped high-Reynolds-number flow, it was 
found necessary to take into account the leading edge 
transitional flow region in order to capture the correct 
physical mechanism that leads to dynamic stall. 
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Introduction 

Application of Navier-Stokes methods to complex 
unsteady, turbulent flows is a field of continuous inter- 
est. An example of such a flow is the massively sepa- 
rated flow over airfoils at high incidence or in dynamic 
motion. For the computation of these flowfields, ap- 
plication of algebraic turbulence models, such as the 
Cebeci-Smith 1 model , the Baldwin-Lomax 2 model, 
the algebraic Renormalization Group (RNG) based 
model 3 or the Johnson-King 4 model, becomes very 
complicated and often ambiguous. The source of am- 
biguity comes from the difficulty in defining charac- 
teristic length scales, e.g., boundary layer thickness 
required by these models. The standard k - <r 5 and 
k — u 6 * * * * * models and the field-equation, (one and two-- 
equation) turbulence models developed over the last 
few years (Refs. 7, 8 and 9), do not have these am- 
biguities. The more recently developed models 7,8,9 
also seem to show promise for massively separated 
flowfields. These models have been extensively tested 
for attached and mildly separated steady flows. The 
ability of these models to resolve massively separated 
flows, such as flow over a wing at stall, a slender body 
at an angle of incidence, or unsteady separated flow 
at dynamic stall conditions, however, has not been 
systematically evaluated. The objective of this inves- 
tigation is to test one and two-equation models for 
unsteady massively separated flows. 

Some of the recently developed turbulence models 
have been tested for two dimensional unsteady sep- 
arated flows. It was found that it is advantageous 
to use the Johnson-King’ model for separated flows , 
where nonequilibrium effects are important. In Refs. 
10-12 it was shown that for light dynamic stall, the 
Johnson-King model captures flow separation bet- 
ter and yields improved predictions of hysteresis ef- 
fects, compared to standard algebraic models. In the 
present investigation, however, it was found that the 
experimental data 13 used for validation of fully tur- 
bulent solutions in Ref. 10 and 1 1 depend also on 



leading edge transition, which has a significant effect 
on the development of the unsteady flowfield. The 
k — e model has been compared to the Baldwin-Lomax 
model in Ref. 14 and 15 for dynamic motion. In a 
previous investigation 16 an extensive evaluation of al- 
gebraic, half and one-equation models for the predic- 
tion of dynamic stall was conducted. A single model 
that could predict accurately the attached unsteady 
flow, the light- and deep-stall regime was not identi- 
fied. 

In Ref. 16 a central difference numerical proce- 
dure validated with experiments in Ref. 10 and 17, 
was used. It was found, however, that unsteady solu- 
tions obtained with this central differencing schemes 
were sensitive to grid spacing in the near wall re- 
gion and to artificial smoothing. Therefore, here the 
upwind-biased numerical scheme described in detail 
in Ref. 18 is used. One and two-equation models are 
tested first for steady separated flow. The same mod- 
els are then applied to compute unsteady flowfields 
over oscillating airfoils. The ability of one and two- 
equation models to predict unsteady attached flow, 
the light and deep stall regime is evaluated. The 
computed results appear to be grid independent and 
no artificial smoothing is used. Finally, it is shown 
that the leading edge transitional flow has a signifi- 
cant effect on the development of the unsteady flow- 
field about oscillating airfoils. The numerical scheme 
and the turbulence models are briefly described in the 
following sections. 

Numerical Implementation 

The thin-layer approximation of the compress- 
ible, Reynolds- averaged, Navier-Stokes equations for 
body-fitted coordinate system, (£, 77 ), is used. These 
equations are as follows 


d* q -f <9^F d n G = Re~ 1 d v S 


here , q is the conservative variable vector, q = 
[p, pu, pv, e] T , F and G are the inviscid flux vectors, 
and S represents the thin-layer approximation of the 
viscous terms in the normal direction. In the above 
equations all geometrical dimensions are normalized 
with the airfoil chord length c; the density p is normal- 
ized with the free-stream density pool the Cartesian 
velocity components, (it, v ), of the physical domain are 
normalized with the free-stream speed of sound a 0 0 ; 
and the pressure p is normalized with poo- 

The following upwind-biased, factorized, itera- 
tive, implicit numerical scheme is used to compute 
the mean flow. 


[I + ht(V>A+ t + A^A ik )] p 
x[/ + Mv;5+ +A 'Br k 

- 

x(^,r-Qy 

— — [(QiJ fc ~ (Fi + i/2'k — Ff-i/2'k) 

(^gjfc + l /2 ~ ^gjfc- 1 / 2 ) 

— Re 1 h T? (5? jt+1 y 2 - 5’fj t _ 1 / 2 )] 


In this equation, = Ar/A£, etc., A ± — ( dF/dQ ), 
etc., are the flux Jacobian matrices, and A, V, and 6 
are the forward, backward and central difference oper- 
ators, respectively. The quantities Gg k+i/ 2 > 

and Si t jfc+ 1/2 are numerical fluxes. 

The inviscid fluxes F and G are evaluated using 
Osher’s 19 upwinding scheme. The numerical fluxes 
for a third-order accurate upwind-biased scheme are 
given by 


^i+l/2, * - ^i+l/2,* + g AF t-l/2,k + 2AF l + l/2,k 
~ 6 [ Z ^^+ 3 / 2 . Jk 2A ^l + l/2,Jk 

= F(Qi,k , + g |AF + (Ci + i ( r;, Qi,k) 

+2AF + (Qi t k, Qi+i,k) 


-g[AF (Qi,k,Qi+i,k) 
+2A F (Qi+i t k, Qi,k) 


Here, F is the first-order accurate numerical flux for 
Osher’s scheme 19 given by 


F i+l/2,k - 2 


[Qi - fl -1 

Fi,k + Fi+i,i- j {P? - F-}dQ\ 


where F, = F+ + F~,F± = (|§) ± , and A F ± are 
the corrections to obtain high-order accuracy. For 
the linearization of the left-hand side terms, the flux 
Jacobian matrices A, B are evaluated by the Steger- 
Warming 20 flux-vector splitting. The viscous fluxes 
■5*, jb-f- 1/2 are computed with central differences. 


Turbulence Models 

An attractive feature of one- and two-equation 
models is that they can be utilized in a more straight- 
forward manner compared to algebraic models in both 
structured and unstructured flow solvers which are 
becoming increasingly more popular. The accuracy 
and the numerical robustness of these models should 
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be^ further demonstrated. Among the most widely 
used one-equation models are the Baldwin-Barth 7 
(B-B) and the Spalart-Allmaras 8 (S-A) models. The 
first model was derived from the two-equation k - e 
model by introducing some simplifying assumptions. 
The Spalart-Allmaras model was developed based on 
dimensional analysis and empirical criteria. An ad- 
vantage of the above field-equation turbulence mod- 
els compared to the algebraic and half-equation mod- 
els is that they do not require evaluation of ambigu- 
ous length scales. One-equation models require nu- 
merical solution of only one partial differential equa- 
tion; therefore, they are less computationally inten- 
sive compared to two-equation models. In this pa- 
per the standard versions of the Baldwin-Barth and 
Spalart-Allmaras models are used to compute steady 
and unsteady separated flows. 

The most popular non-algebraic turbulence mod- 
els are two-equation eddy-viscosity models. These 
models solve two transport equations, one for the tur- 
bulent kinetic energy k and another one related to the 
turbulent length- (or time-) scale. Among the two- 
equation models, the k — e model is the most widely 
used today. The original Jones and Launder k — e 
model 5 and its variations have been very successful 
in a large variety of different flow situations, but it 
also has a number of well known shortcomings. From 
the standpoint of aerodynamics, the most disturbing 
is the lack of sensitivity to adverse pressure-gradients. 
Another shortcoming of the k — e model is associated 
with the numerical stiffness of the equations when in- 
tegrated through the viscous sublayer. 

The k — co model has been developed by Wilcox 6 
to overcome the shortcomings of the k — e model. This 
model solves one equation for the turbulent kinetic 
energy k and a second equation for the specific tur- 
bulent dissipation rate (or turbulence frequency) u ;. 
The k — co model performs significantly better under 
mild adverse pressure-gradient conditions than the 
k — e model. Another strong point of the model is 
the simplicity of its formulation in the viscous sub- 
layer. The model does not employ damping functions 
and has straightforward boundary conditions. This 
leads to significant advantages in numerical stability. 
The k — co model has been validated extensively 6,21 
for many flow cases with and without adverse pres- 
sure gradient. For all cases it was found to perform 
equally well or better than the k — e model. The major 
shortcoming of the k — co model is that the results of 
the model depend strongly on the freestream values, 
coj, that are specified outside the shear-layer. 

The free stream dependency of the original Wilcox 
k — to model has. been investigated in detail in Ref. 
22, and it has been shown that the magnitude of the 


eddy viscosity can be changed by more than 1 00% 
just by using different values for co f . This is clearly 
unacceptable and corrections are necessary to ensure 
unambiguous solutions. The standard k - to model 
developed by Wilcox has been modified in Ref. 9 
so that the computed solutions are insensitive to the 
freestream values oft oj . This modified model is called 
Baseline (BSL) k— co model. The BSL k— co model was 
further modified 9 in order to improve the predictions 
of strong adverse pressure gradient separated flows, 
this model is called Shear Stress Transport (SST) 
k — co model. In this paper, the SST k — co turbu- 
lence model will be extensively tested for unsteady 
flows. This model has been tested in Ref. 9 for a 
wide class of steady separated flows and has shown 
good agreement with experiments. 


Baldwin-Barth (B-B) model 


The eddy viscosity of the Baldwin-Barth one- 
equation model 7 is given by u t = uc^f^Rr = vc^^Rt, 
here R? is the turbulent Reynolds number and R t the 
modified turbulent Reynolds number. The quantity 
Rt is the solution of the following field equation 


=( c ej2 ~ c €l )\fvR T P 

-f- (z/ 4- — )V 2 (i^^7’) 

— — (Vzq) * V(uR t ) 


This is a partial differential equation for the field 
quantity R T = k 2 /ue = Rrh{Rr ), and 


— =( c < 2 - c fl ) v ^/ K2 
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/ 2 (y + ) 
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[s/^dC + ^==(Aeip(-y+/>l+)D 2 
+ -^+exp(-y + /At)Di)} 



where y + — u T y/u and u T is the skin friction velocity. 
The constants of the model are: 

k =0.41, c fl = 1.2, c i2 — 2.0 

Cf t =0.09, A + = 26., A 2 + = 10. 

This model is applied to the entire flowfield to com- 
pute the eddy viscosity. 


Spalart-Allmaras (S-A) model 

The second one-equation model used is the Spalart- 
Allmaras (S-A) model 8 . The eddy viscosity is ob- 
tained from the solution of the following partial dif- 
ferential equation for u. 


DD 

~Dt 


=Q>i(l - fn)Sv 


1 l 

+ — V • ((*/ + i>)Vv) + cj 2 (Vv) 2 


f C &1 , 

Cwlju) o Jt2 


IT 2 


d\ + ftlAU2 


here S is the vorticity magnitude and S = S + 

7TT f v 2 1 /v 2 ~ 1 — l+x/ui ’ f vl = 1 " l + y3+ Cul . c v\ = 

7.1 and d is the distance to the closest wall. The other 
functions of the model are: 

fti =c n g t exp[ - c t 2 ^jji[d 2 + g'fd 2 ) J 

ft2 -c iZ exp(-c t 4X 2 ) 

g t =mm(0.1, AU/u>tr Az) 


where x = v/v and u> tr is used here to denote the 
vorticity at the wall at the boundary layer trip point. 
The constants of this model have been chosen the 
same as in the original reference 8 . The constants of 
the model are: 


c bl =0.1355, c 62 = 0.622, <r = 2/3 

c bl , 1 + Cj, 2 


Cwi ~~ + 


Cw2 — 2.0, C-w3 — 0.3 


9 + 

V\ 


ff =r + o„ 2 (r -r),r= SkH , 

Cil =1.0, Ct2 = 2.0, Ct3 = 1.1, C<4 = 2.0 


The eddy viscosity is given as 


”t - vfvi 

and the Reynolds shear stress are given by —upuj = 
2 /'t Ab- 


original k — u) model 


The original k — cj model is given by 


Dpk 

~dT 

Dpu> 


= r, 


dui 
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” axj ~ + 


P + Pt) 


dk 

dxj 
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w 




dui 2 d f / iv \ doj 


Dt v t 

The constants of the original Wilcox model are 
<rf =0.5, cr^ = 0.5, 0 W = 0.0750, 

0" =0.09, k = o.4i, J w = p w //?* 


BSL fc — model 

The (BSL) model is identical to the cj model of 
Wilcox 6 for the inner region of a boundary layer (up 
to approximately 8/ 2) and gradually changes to the 
standard k — € model in the outer wake region. In or- 
der to enable computations with one set of equations, 
the k — € model was first transformed into a k — oj 
formulation. The blending between the two regions 
is performed by a blending function that gradually 
changes from one to zero in the desired region. No a 
priori knowledge of the flowfield is necessary to per- 
form the blending. The function also ensures that the 
k — e formulation is selected for free shear layers. The 
performance of the new (BSL) model is very similar to 
that of the Wilcox k — u> model for adverse pressure 
gradient boundary-layer flows (and therefore better 
than that of the k — e model), but without the unde- 
sirable freestream dependency. For free shear layers 
the new model is basically identical to the k — e model, 
which predicts spreading rates more accurately than 
the k — u model. The Baseline (BSL) k — u> model is 


Dpk 

~W 

Dpoj 

~Jn 


•ij 


dxj 


■ p pUJK 


dxj 

d 






dxj 

du> 


7 du{ 2 d N du 

u) dxj dxj 


where the constants of the model are computed as <p = 
Fi<j) W 4- (1 — F{)4> K where F\ is a blending function 
as defined in Ref. 9 and <j) W and <f> K the constants for 
the original k — u and the k — e model, respectively. 
The following standard values for the k — e model are 
used, 


<rf =1.0, <t% = 0.856, 0 K = 0.0828, 

/?• =0.09, k = 0.41, j K = /?*7/r - 

corresponding to the constants C t \ = 1.44, C c2 = 1.92 
of the k — e model. 
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k — € model 

The fc — € model implemented here is based on 
the same formulation as the BSL model, except that 
the switch from the Wilcox model, constants <j> w , to 
the i k — e model, constants <fr K , takes place not in the 
wake region of the boundary layer but just outside the 
viscous sublayer. With this formulation, the Wilcox 
model is only used as a sublayer model and the model 
is referred to as a two-layer k — € model. 

SST k — u> model 

Although the original and the new BSL k — u> 
model perform better in adverse pressure gradient 
flows than the k - e model, they still underpredict 
the amount of separation for severe adverse pressure 
gradient flows . 9 In an attempt to improve matters, 
the eddy-viscosity formulation of the BSL model is 
modified to account for the transport effects of the 
principal turbulent shear stress. The motivation for 
this modification comes from the Johnson-King (J~ 
K) model 4 which has proven to be highly successful 
for adverse pressure gradient flows. The J-K model 
is based on the assumption that the turbulent shear 
stress is proportional to the turbulent kinetic energy 
in the logarithmic and wake regions of a turbulent 
boundary layer. Johnson and King solve an equa- 
tion for the maximum turbulent shear stress at each 
downstream station and limit the eddy viscosity in or- 
der to satisfy this proportionality. In the framework 
of two-equation models, the turbulent kinetic energy 
is already known and it is therefore only necessary to 
limit the eddy viscosity to account for the same effect. 
The resulting model is called Shear Stress Transport 
(SST) model. For the SST k-w model, the constants 
<f) W of the BSL model are replaced by the constants 
< fi s as follows: 

4 =0.85, 4 = 0.5, f = 0.0750, a = 0.31 

/T =0.09, k = 0.41, 7 s = 4 IP* - 4« 2 /4F 

where the same convention 4> = Fi<f> s + (1 — F\)<j> K is 
used and the eddy viscosity is given by 

a\k 

Vt max(au>] QF 2 ) 

Where O is the vorticity magnitude and F 2 = tank 
{argl) with arg 2 = ^£ 7 ) 


Results and Discussion 
The main objective of the present investigation 
is to assess the accuracy of one- and two-equation 


turbulence models for the computation of unsteady, 
massively separated, high Reynolds number airfoil 
flows. Experimental measurements for dynamic stall 
are usually available in the form of unsteady loads 
(lift coefficient Cj, drag coefficient C d and pitching 
moment coefficient C m ). Accurate computation of 
these quantities for unsteady flow is important, and 
usually is sufficient for engineering and design appli- 
cations. Sometimes unsteady surface pressures are 
also reported. 

To confirm that the turbulence models under con- 
sideration had been implemented correctly, computa- 
tions were first performed for a standard incompress- 
ible airfoil test case 23 with trailing edge separation. 
Solutions over the NACA-4412 airfoil are obtained at 
a = 13.87° angle of attack and Re c = 1.52 x 10 6 . In 
the experiment of Ref. 23, the measurements are ob- 
tained for low speed incompressible flow. Here, the so- 
lutions with different turbulence models are obtained 
at Moo — 0-2. The computed surface pressure coef- 
ficients are compared with the measurements of Ref. 
23 in Fig. 1. 



x/c 


Fig. 1 Comparison of the computed and measured 
surface pressure coefficient with the experiment of 
Ref. 24; a = 13.87°, Re - 1.52 x 10 6 . 

Only small differences in the computed surface pres- 
sure coefficients are obtained in the trailing edge 
region. The computed velocity profiles at various 
streamwise locations are compared with the measured 
values in Fig. 2. 
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Fig. 2 Comparison of the computed and measured 
boundary layer profiles with the experiment of Ref. 
24 for x/c = 0.675, 0.731, 0.786, 0.842, 0.897, and 
0.935; a = 13.87°, Re = 1.52 x 10 6 . 


The B-B model overpredicts separation and the S-A 
model yields less separation compared to the experi- 
ment. The k—e and the k— u models predict attached 
flow. The SST k — co model predicts the separation 
and profile shape fairly well. The trends observed for 
this separated flow at a fixed angle of attack carry 
over to the unsteady case and help to interpret the 
computed unsteady results. 


Fully Turbulent Tripped Unsteady Flows 

For validation of unsteady, fully turbulent solu- 
tions, the experimental measurements of Ref. 24 for 
a NACA-0015 airfoil are used, because in this exper- 
iment, as opposed to that of Ref. 13, the boundary 
layer was tripped at the leading edge to ensure a fully 
turbulent boundary layer for attached, light stall and 
deep stall cases. Oscillatory motions at sufficiently 
high angles of incidence include both massive separa- 
tion during the upstroke and flow reattachment dur- 
ing the downstroke. The free stream Mach number is 
Moo =0.3 and the Reynolds number, based on the 
airfoil chord length is, Re c = 2 x 10 6 . In the exper- 
iment the flow was tripped at the leading edge and 
it is expected that the surface flow is fully turbulent. 
The airfoil oscillates as a(t) — a m + a a sin(cot) with 



a reduced frequency k = 0.1. The oscillation ampli- 
tude remains fixed at o a = 4.2° and variation of the 
mean angle a m leads to different flow regimes. At- 
tached flow corresponds to a m =4°. The light and 
deep stall regimes are obtained for a m = 11°, and 
a m — 15°, respectively. 

The computations use a C-type 311 x 91 point 
grid with 130 points on the suction side and 45 points 
in the wake having a grid spacing dz = 0.00001 chord 
lengths away from the airfoil surface as baseline grid. 
Normal spacing of 0.00001c yields a y + ss 2 for the 
first grid point above the suction side surface. A 
421 x 151 points grid with increased resolution in the 
normal to the wall separated flow region and reduced 
normal grid spacing of dz — 0.000005, is also used. 
This grid has 111 points on the pressure side and 221 
points on the suction side. Two oscillatory cycles are 
computed for every case and the third cycle is always 
identical to the second cycle. Different number of 
time steps for a cycle was tested, and it was found that 
computations with 10000, 16000 and 40000 time steps 
gave identical solutions. The results presented are ob- 
tained using 16000 time steps, which correspond to a 
nondimensional time step At « 0.0065 or a Courant 
number Cu « 700. The performance of all turbulence 
models is evaluated for the deep stall case. For the 
light stall case, solutions are computed only with the 
B-B, S-A and SST k — co turbulence models, and a 
grid resolution study is performed. The attached un- 
steady flow case it is computed only with the B-B 
and SST k — co models only. 

a(t) = 4° +4.2 °sin(t) 

This flow is essentially attached and it is com- 
puted only with the B-B and SST k — co models. 
Comparisons of the computed hysteresis loops with 
the experiment (Fig. 3) shows that the loads com- 
puted with the SST k — co model are in good agree- 
ment with the experiment. The lift computed with 
the B-B model overpredicts the experimental values 
and the pitching moment hysteresis does not agree 
with the experiment. The computed drag, however, 
is in agreement with the experiment. Comparison of 
the unsteady surface pressure coefficients (Fig. 4) at 
a — 4° and a = 8° during the upstroke computed 
by the two models, shows that the lift and pitching 
moment disagreements are not caused by large dif- 
ferences of the computed unsteady surface pressures. 
Due to the small trailing edge separation, the varia- 
tion of the pitching moment and drag coefficients is 
very small (an order of magnitude less than that ob- 
served in the deep stall case). 
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Fig. 4 Comparison of the surface pressure coefficient 
at a = 4° up and a = 8° up computed with the B-B 
and the SST k — u> models. 


a(t) = 11° + 4.2 °sin(t) 

The computed solutions show that this flow is 
characterized by moderate trailing edge separation 
which develops at the peak of the cycle. The flow 
remains separated for a large portion of the down- 
stroke and a recirculatory region of about half a chord 
length is observed. As a result, more significant hys- 
teresis effects than the previous case are obtained. 
The hysteresis loops obtained from solutions using the 
baseline grid with the B-B, the S-A, and SST k — w 
turbulence models are compared with the experimen- 
tal data in Fig. 5. The B-B model predicted the 
most separation and yields a lower lift during reat- 
tachment but it gives good predictions for drag and 
pitching moment coefficients. The S-A model pre- 
dicts the least separation and shows earlier flow reat- 
tachment. Similar to the S-A model, the SST k — w, 
model even though it shows closer agreement with the 
experimental lift during the initial part of the down- 
stroke, predicted more rapid lift recovery. A solution 
computed on the refined 421 x 151 point grid with 
the SST k — ui turbulence model is compared with 
the baseline grid solution in Fig. 6. Little grid sensi- 
tivity is obtained with the refined grid. 

a(t) = 15° + 4.2 °sin(t) 

The computed solution show that this unsteady 
flowfield is characterized by massive flow separation 
which develops before the peak angle of incidence. 


At peak incidence and before the downstroke the dy- 
namic stall vortex is shed and a trailing edge vortex 
forms. Shedding of the dynamic stall vortex causes 
decrease in lift and pitching moment. The flow re- 
mains separated for a large part of the downstroke 
and significant hysteresis effects are obtained. It is 
a challenge for turbulence models to be able to cap- 
ture flow separation and reattachment and yield good 
quantitative predictions for these kind of flows. The 
lift, drag, and pitching moment hysteresis loops ob- 
tained from computations with the B-B and S-A one- 
equation turbulence models are compared with the 
experiment in Fig. 7. The lift hysteresis is predicted 
reasonably well by both models. The drag and pitch- 
ing moment hysteresis loops indicate that both mod- 
els delay onset of separation. The loads computed 
with the B-B model show oscillatory behavior at the 
downstroke. The B-B model predictions are in closer 
agreement with the experiment. A smaller extent of 
separated flow is obtained with the S-A model, re- 
sulting in smaller extreme values of drag and pitch- 
ing moment. A more rapid flow reattachment was 
also observed. 

Predictions of hysteresis loops obtained from so- 
lutions with the two-equation turbulence models are 
compared with the experiment in Fig. 8. The k - e, 
the original k — w, and the BSL k—u (where the free 
stream dependency is removed) did not yield enough 
separation. Therefore, the loads computed with these 
models significantly deviated from the experimentally 
measured values. The solution obtained with the SST 
k — cj model shows large flow separation and the pre- 
dictions are in close agreement with the experiment. 
However, at large angles of incidence and during the 
downstroke the loads show again oscillatory behavior. 

Solutions computed with the B-B and SST k — co 
models show development of a trailing edge vortex at 
the end of the pitch-up motion. During the down- 
stroke this vortex convects in the wake and another 
trailing edge vortex forms. The second vortex ini- 
tially grows in size and then convects in the wake and 
the vortex shedding repeats. This vortical activity 
and the suction side separated flowfield is shown by 
a series of snapshots during the downstroke in Fig. 
9. In Fig. 9, the vorticity magnitude computed us- 
ing the B-B model is used to describe the flowfield. 
During the downstroke the trailing edge vortex shed- 
ding causes oscillations to the loads. The computed 
solutions for this part of the cycle are sensitive to grid 
distribution at the trailing edge and the wake region. 

The unsteady solutions follow the trends of the 
steady results. The B-B model has the tendency to 
predict the most separation following the SST k - u 
model. The S-A model tends to underpredict the 
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amount of separation but not nearly to the degree of 
the standard k — u> and k - e models. To better put 
the results in perspective, the pitching moment for all 
three oscillatory cases are plotted to the same scale 
in Fig. 10 and compared to the B-B and SST k — u 
model solutions. The drag and pitching moment vari- 
ation for the deep stall case is an order of magnitude 
larger compared to the variation in the attached flow 
and the light stall flow cases. 
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frequency k = 0.1, M <*> = 0.3 with untripped flow of a 
Reynolds number Re c — 4.0 x 10 6 . The same case was 
also considered in Refs. 10 and 11. However, because 
significant hysteresis effects were not obtained for the 
experimental maximum angle of incidence, the ampli- 
tude of the oscillation or the mean angle were slightly 
increased for the computation. As a result, in both 
investigations a maximum angle of attack larger than 
a ma x = 15° reported by the experiment was reached. 
Experimental uncertainties and tunnel wall interfer- 
ence effects justified this alteration of the experimen- 
tal conditions. Another reason why the maximum 
angle of incidence had been increased was to promote 
separation predicted by the turbulence models. Hys- 
teresis effects were obtained at these larger angles of 
incidence, and it was concluded that massive flow sep- 
aration at the trailing edge alone was responsible for 
stall. 


o.i 
o.o 
o -o.i 
- 0 . 2 - 
- 0.3 

0 5 10 15 20 

Angle of Attck, deg. 

Fig. 10 Effect of increased mean angle on pitching 
moment. 

As seen in Fig. 10, the discrepancies between experi- 
ment and computation for the attached and light stall 
cases are very small in relation to those for the deep 
stall case. 

The effect of transition 

A light stall case from the experimental results re- 
ported in Ref. 13 is chosen to demonstrate the effect 
that leading edge transition can have on the develop- 
ment of the unsteady flowfield. In Ref. 13, measure- 
ments in the form of integrated aerodynamic loads 
( Ci,Cd,C m ) have been reported for a wide range of 
flow conditions and airfoil shapes. In addition, un- 
steady surface pressure coefficients are given. .In con- 
trast to the NACA-0015 experiment of Ref. 24, in 
Ref. 13 boundary layer trips were used only for lim- 
ited number of deep stall cases. Unsteady solutions 
are obtained for a NACA-0012 airfoil executing har- 
monic motion a(t) — 10° 4- b°sin(ujt), with a reduced 



Fully turbulent flow simulations of the previous 
section have demonstrated that the B-B model pro- 
duces the most separation. The S-A model produces 
less separation compared to the B-B and the SST 
k — uj models. The B-B and the SST k — to mod- 
els yield similar predictions, but the B-B model is 
more computationally efficient. Therefore, this case is 
solved only with one-equation models. A solution ob- 
tained with the B-B model, which produces the most 
separation, for the same oscillation amplitude as the 
experiment 13 still did not yield significant hysteresis 
effects. In addition, a counterclockwise loop for the 
pitching moment, as opposed to the clockwise loop 
shown in the experiment, was obtained: The same 
procedure of Refs. 10 and 11 was followed again, and 
the oscillation amplitude was “arbitrarily” increased. 
A solution with the B-B model was obtained with 
an oscillation amplitude of 5.3°. For the S-A model, 
which yields less separation, the oscillation amplitude 
was further increased to 5.5°. The loads obtained 
with the B-B and the S-A models are compared with 
the experiment in Fig. 11. The lift and pitching mo- 
ment hysteresis reasonably agree with the experiment. 
The computed drag, even though it follows the exper- 
imental trends, underpredicts the extreme measured 
drag values. Also the phase angle where the computed 
drag and nose down pitching moment increases, lags 
the experimental values by approximately one degree. 
Comparison of the measured unsteady surface pres- 
sure measurements with the computed surface pres- 
sure, indicates that large discrepancies occur on the 
suction side, and that the agreement of the lift and 
the pitching moment with the experiment is coinci- 
dental. 

Careful observation of the experimental surface 
pressure measurements shows that the lift drop, the 


8 



increases in drag, and nose-down pitching moment 
is associated with a drop in the leading edge suction 
peak. This suction pressure drop occurs because of 
leading edge flow separation. The two physical mech- 
anisms that can force the flow to separate at the lead- 
ing edge are either a shock or laminar/transitional 
flow behavior. The progressive drop of the suction 
peak shown in the experiment does not support the 
shock separated flow assumption. Therefore, there is 
an indication that, even though the Reynolds number 
is large, the leading edge flow could be transitional. A 
rough approximation of the transitional flow behav- 
ior at the leading edge is performed with the following 
procedure. The transition onset is specified to occur 
immediately downstream of the suction peak location. 
The flow from the stagnation point until the transi- 
tion onset is computed els laminar. The production 
term of the one-equation B-B model is set equal to 
zero for the laminar region. As a result, the model 
yields an eddy viscosity of almost zero for the laminar 
region. Beyond the transition point, the full produc- 
tion term is utilized and the computed eddy viscos- 
ity rapidly increases downstream from the transition 
point until a fully turbulent value is reached. 

A 351 x 91 point grid with 170 points on the suc- 
tion side is used for the numerical solution. This 
grid has refined resolution at the leading edge re- 
gion. In the computed solutions the transitional 
flow region extends only over a few streamwise com- 
putational cells. A leading edge separation bubble 
forms at approximately 14 degrees during the up- 
stroke and increases significantly in size before the - 
peak of the cycle. The loads obtained from the lam- 
inar/transitional/turbulent flow solutions are com- 
pared with the experiment in Fig. 12. For compari- 
son, the loads obtained from a fully turbulent solution 
are shown in the same figure. The lift hysteresis loop 
(Fig. 12 a) shows good qualitative agreement with 
the experiment. Good quantitative agreement is ob- 
tained for the upstroke but a more rapid lift recovery 
during part of the downstroke is observed. Similar 
trends are shown for the drag and pitching moment. 
It appears that the transitional solution predicts a 
more rapid flow reattachment. However, there is good 
agreement with the experiment for the computed nose 
down pitching moment and drag increase attributed 
to the massively separated flow during the initial part 
of the downstroke. It is also significant that the ex- 
treme values of the drag and pitching moment are 
closely predicted and the computed loads do not lag 
the experiment. Discrepancies obtained for the down- 
stroke are caused by uncertainties in transition mod- 
eling and deficiencies of the turbulence model. 


The surface pressure coefficient distributions for 
three angles during the upstroke, a = 14.0°, a = 14.5° 
and a = 14.9°, obtained from the fully turbulent and 
the transitional computations are compared with the 
measured values in Fig. 13. At a = 14.0° (Fig. 13 
a) both the fully turbulent and the transitional so- 
lutions are in agreement with the experiment for the 
region near the leading edge. The fully turbulent solu- 
tion, however, slightly overpredicts the suction peak. 
The surface pressure distribution obtained from the 
transitional solution is in close agreement with the 
experiment and shows the development of a leading 
edge separation bubble. A very small region of tran- 
sitional flow is found. As the angle of attack increases 
to a = 14.5°, the experiment shows a large drop in 
the suction peak caused by flow separation. Visu- 
alization of the computed velocity fields shows that 
only trailing edge separation was obtained for the 
fully turbulent computation. The transitional solu- 
tion, on the other hand, yields, more separated flow 
and shows formation of a vortex-like structure. At 
a = 14.9°, which is the peak angle of attack of the 
experiment, the suction peak is further diminished 
and the vortical structure generated at the leading 
edge is convected downstream. At this angle the fully 
turbulent solution does not agree with the experiment 
and shows further increase of the suction peak. It ap- 
pears that the transitional solution properly captures 
the physical mechanism observed in the experiment 
and shows development of a leading edge, vortical, 
dynamic-stall-like structure as the computed surface 
pressure coefficient demonstrates. 


The values of lift and pitching moment obtained 
from the turbulent solution with an increased oscil- 
lation amplitude (Fig. 11) are coincidentally close to 
the measured values. The computed drag coefficient, 
however, disagrees with the experiment. Therefore, it 
is necessary to always compare lift, drag and pitching 
moment coefficients, when surface pressure measure- 
ments are not available. The fully turbulent solution 
predicts attached leading edge flow and trailing edge 
separation. The transitional solution predicts a lead- 
ing edge vortex-like structure and larger overall sepa- 
ration, compared to the fully turbulent solution. The 
leading edge flow development affects significantly the 
suction surface fiowfield and results in larger overall 
separation. The leading edge vortical structure is not 
the same as a classical dynamic stall vortex, which 
is clearly observed in both experiments and compu- 
tations for larger oscillation amplitudes or different 
pitch rates. This structure forms at about 14 degrees 
during the upstroke and bursts in the boundary layer. 
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Around the peak angie of the cycle a rapid progres- 
sion of the trailing edge separation towards the lead- 
ing edge is also observed. 

As an approximate location for the transition on- 
set, the maximum suction pressure point is used. For 
high Reynolds number flow, the transitional region is 
small. It appears that in the computation, where an 
approximate transition onset location is determined 
based on the maximum suction pressure, and a sim- 
ple transition model that yields an “effective” eddy 
viscosity for the transitional region is used, do not sig- 
nificantly degrade the solution for the pitch-up part 
of the cycle. It is not expected that these rough ap- 
proximations allow accurate modeling of the complex 
physical mechanisms of transition, bubble formation 
and reattachment. It is demonstrated, however, that 
transition plays a significant role in the development 
of the unsteady separated flowfield. 

Conclusions 

An evaluation of the ability of one- and two- 
equation turbulence models in predicting hysteresis 
effects of unsteady fully turbulent flow over oscillat- 
ing airfoils in the light and deep stall regime was con- 
ducted. None of the models considered in this inves- 
tigation is capable of accurately predicting the deep 
stall case. However, the B-B, the S-A, and the SST 
k — to models show a significant improvement over 
standard two-equation models. For the light stall 
case the S-A model did not yield sufficient separa- 
tion and underpredicted the extreme values of the 
unsteady loads. The B-B model overpredicted the 
lift drop for the light stall case and the attached flow 
cases. The standard k — e and the k — to models did 
not predict separation even for the deep stall case. 
The SST k — u model gave good predictions for the 
attached and the light stall cases. 

It was found that the leading edge transitional 
flow is of primary importance to the overall devel- 
opment of the unsteady flowfield, if the flow is not 
tripped at the leading edge. A laminar/transitional 
leading-edge separation bubble developing during the 
pitch-up motion produces a dynamic-stall-like vorti- 
cal structure. It was shown that a simple transition 
model significantly improves the predictions. How- 
ever, accurate methods for transition modeling and 
prediction are still required. 
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Fig. 5 Hysteresis effects for light stall flow; M — 0.3, a(t) = 11° + 4:2°sin(t), k - 0.1, 
Re = 2 x 10 6 (fully turbulent). 
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Fig. 6 Effect of grid refinement on the computed hysteresis effects; 

M ~ 0.3 ,(x(t) = 11° + 4.2°sm(i), k = 0.1, Re — 2 x 10 6 (fully turbulent) 
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■ r Hysteresis effects for deep stall obtained with one- equation models; 

M = 0.3, a(t) = 15° + 4.2 °sin(t),k == 0.1, Re = 2 x 10 6 (fully turbulent) 






0 - 1-1 


o.oi 


E 

O 


- 0.1 1 


-0.2 H 


(c) 



- 0.3 


10 


“I — 
12 


14 16 • 18 

Angle of Attck, deg. 


20 


Fig. 8 Hysteresis effects for deep stall obtained with two-equation models; 

M = 0.3, Oi(t) = 15° +A.2°sin(t), k = 0.1, Re = 2 x 10 6 (fully turbulent). 
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Fig. 11 Effect of increased oscillation amplitude; M = 0.3, a(t) — 10° + 5°sm(t), k 
Re = 4 X 10 6 (fully turbulent). 




O Experiment, McCroskey 

Computed, B-B 

Computed, B-B 

model 

mode! 








Fig. 13 Effect of transition on the computed surface pressure coefficient; 
M = 0.3, oc(t) = 10° + 5 °sm(i), fc — 0.1, Re = A x 10 6 
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